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Qh We report on an improvement to the implementation of the Maximum Entropy Method (MEM). 

Q_i It amounts to departing from the search space obtained through a singular value decomposition 

(SVD) of the Kernel. Based on the shape of the SVD basis functions we argue that the MEM 
spectrum for given N T data-points D(t) and prior information m(co) does not in general lie in 
this N x dimensional singular subspace. Systematically extending the search basis will eventually 
recover the full search space and the correct extremum. We illustrate this idea through a mock data 
analysis inspired by actual lattice spectra, to show where our improvement becomes essential for 



the success of the MEM. To remedy the shortcomings of Bryan's SVD prescription we propose to 
i use the real Fourier basis, which consists of trigonometric functions. Not only does our approach 

lead to more stable numerical behavior, as the SVD is not required for the determination of the 
basis functions, but also the resolution of the MEM becomes independent from the position of the 
reconstructed peaks. 
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1. Introduction 

Lattice QCD is a central pillar for the non-perturbative investigation of the strong force. One 
important aspect is that its formulation in Euclidean time is amenable to Monte Carlo simulations 
at arbitrary T , i.e. also close to a phase transition. Computational advances over the last decade 
have furthermore made it possible to consistently probe the properties of e.g. the high temperature 
state of matter, the quark-gluon plasma, with staggered dynamical quarks in the continuum limit 
[1] and to deploy the even more costly fermion formulations, such as Wilson or domain wall [2]. 

One particular area, which directly benefits from the availability of non-perturbative infor- 
mation about the QCD medium is the study of heavy quarkonia and its proposed melting [3] in 
relativistic heavy ion collisions. In these studies the dynamics of a pair of a heavy quark and anti- 
quark is formulated in an effective field-theory framework (for a brief overview see e.g. [4] and 
references therein), which allows to make a close connection to the physics of the medium by re- 
lating the non-relativistic static potential governing the dynamics of the heavy QQ to the real-time 
rectangular Wilson loop Wu (r,t) = (exp — if n dx^A^ \ . 

Unfortunately, formulating lattice QCD in imaginary time deprives us from direct access to dy- 
namical information such as real-time observables. Field theory however provides a handy concept 
to relate imaginary and real-time quantities, through the use of spectral functions. These positive 
definite quantities p(ft)) are connected to the observables D(z) via an integral transform 

/oo 
K(z,co)p(co)dco (1.1) 
-oo 

with analytically known Kernel K(z,co). Since all time dependence is explicitly located in the 
function K(z, to) the analytic continuation z — > it is straight forward, once the values of p(ft)) are 
determined. (In the following we will focus on the exponentially damped Kernel K(z, ft)) = e~ mt ). 

In practice one obtains from simulations of lattice QCD a noisy estimate of the observable 
D(z) at a discrete number of A T datapoints D,. Extracting a function p/ at N m S> N T frequencies 
from such a limited data-set via the inversion of 



D/ = £*bP/Aa* (1.2) 
i=i 



is an ill-defined problem. How to give meaning to such a task is the central question elucidated 
by the Maximum Entropy Method [5, 6, 7, 8]. This branch of Bayesian inference generalizes and 
ameliorates the concept of % 2 fitting by emphasizing the use of so called prior information. I.e. one 
introduces a function m(ffl) which per definition is the correct spectral function in the absence of 
measured data 1 . 

From Bayes theorem we learn that the probability of a test spectral function to be the correct 
spectral function P[D\p]P[p\m] 

P M EM[p\D,m}= p[^| m ] ( L3 ) 

'A complication exists in that since our prior information is usually derived from perturbation theory at high en- 
ergies, there is essentially no information available on the low CO part of the spectrum, we are interested in. Thus we 
already know that the prior function used is not the correct solution in the absence of data. 
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for given N T measured data D(t,-) = A and N w prior information m{(Oi) = mi can essentially be 
written as a product of the likelihood probability 

(1.4) 



P[D\p]=e-^ = ^ V \-\ £ (Di-^Cr/iDj-D?] 
with Cjj denoting the covariance matrix and the prior probability 



r Nm ( Pi 

P[p\m] =e ay = exp aY (pi-m,- p, log [— 



(1.5) 



The parameter a introduced in eq.(1.5) is selfconsistently determined and does not enter the final 
result [6, 7]. The task to perform is to find the global extremum of eq.(1.3), 



g-ftmt\p\D,m] 
bp 



= (1.6) 

P=Pmem 



which is in general very costly, since it involves an Na> dimensional search space. Note that it has 
been proven in [8] that the solution of eq.(1.6) is unique if it exists. This statement does neither 
depend on a particular parametrization of the test spectral function p(co) nor on the number of 
available data-points. Indeed in the extreme case of zero data-points the solution spectrum of the 
MEM would be the prior function. 

In maximizing Pmem[p\D,ii] the likelihood and the prior term compete. Pure % 2 fitting is 
under-determined, so an infinite number of solutions exist, which all reproduce the measured data 
within the errorbars. The prior then selects from within these solutions the one that is closest to 
m((o) while still respecting the constraints from the data. In essence one regularizes the % 2 fitting, 
with parts of the spectrum being determined by the data, parts of it being fixed by the prior. Thus 
to find out which part of the spectrum is actually encoded in the lattice QCD correlator, we need to 
redo the MEM with several distinct prior function to observe which part of p MEM ((o) is invariant 
under a change of m(G>). 

To carry out the maximization procedure in eq.(1.6), the state of the art implementation of the 
MEM relies on the singular value decomposition prescription introduced by Bryan [5], (For an 
alternative implementation without the use of the SVD see [9]) In the following we will argue that 
the global extremum of eq.(1.3) does not in general lie in his N T dimensional singular search space 
but one needs to systematically approach the full N a dimensional space. 

2. Extended Search Space 

Let us quickly revisit Bryan's argument by inserting eq.(1.4) and eq.(1.5) into eq.(1.6) and 
making explicit the positive definiteness of the spectrum by using p/ = m/exp[a/], which together 
with the SVD of K< = VLV 1 yields 

-aa = K t ^- = UW t ^-. (2.1) 
dDP dDP 

Since E is a N m x N m diagonal matrix with only A^ T entries different from zero, Bryan concluded that 
eq.(2.1) restricts the values of a to the space spanned by the first N t columns of U. We check the 
validity of this claim by numerically investigating the behavior of the corresponding basis functions 
U'((Oi) = Ua (/ = !,... ,N r ) of Bryan's search space. 
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discretize the interval co G [—10,20], while z £ (0,6.1] is discretized in N T = 12 steps. Note that at posi- 
tive frequencies above CO — 5 no structure remains, so that reconstruction of peaks becomes very difficult, 
(right) The first 50 columns of the matrix U as comparison, where the oscillations are visible over the whole 
frequency range. 



In Fig.l we show the basis functions of the Kernel K(co, z) = e~ ar , where we have chosen to 
discretize eq.(l.l) over a finite interval CO G [dVn = —10, GWx = 20]. The functions U'(coi) exhibit 
a common feature, i.e. they posses an oscillating part, which starts from cOmin and extends up to 
a certain co osz . Beyond this frequency the functions rapidly damp to zero. Hence in the example 
shown, with the choice 0) m ; n = — 10, the first N x = 12 basis function do not allow the reconstruction 
of peak structures in the region CO > co osz ~ 0. The conceptual problem of Bryan's SVD basis lies 
in the fact that the choice of co^m is completely arbitrary. As long as the resolution of the frequency 
interval stays the same Aft) = (<w A ^ a)mm > extending the frequency range to smaller and smaller CDmin 
is admissible. Decreasing G) m i n does not change the length of the oscillatory part it just shifts 
the basis functions to smaller and smaller frequencies. At some point, the MEM in the singular 
subspace cannot reproduce the measured datasets, since it does not have any peak structures left at 
positive frequencies, while in the full search space the data can be reconstructed without problem. 

In other words, no matter how many datapoints have been measured it is always possible to find 
a CDmin so that e.g. the positive frequency range cannot be resolved. We conclude that the success 
of Bryan's MEM prescription depends on the arbitrary choice of co m i„. The proof for existence and 
uniqueness of the MEM solution on the other hand does not require a certain parametrization of 
the spectral function and thus is valid for the full N m dimensional space. This tells us that the SVD 
space artificially restricts the optimization without guaranteeing the presence of the correct global 
extremum. 

Thus as a first step we propose to systematically enlarge Bryan's search space [10] to approach 
the full M. Nm , which has to contain the correct solution. To this end we parametrize the spectral 
function with N ext > N r parameters bj 

-Next 



Pi = m/exp 



I Uijbj 

7=1 



(2.2) 



As a rule of thumb we can stop increasing N ext once the value of PMEM\p\D,m] stops to grow after 
adding another basis vector. 

To illustrate how the extension of the search space can improve the MEM reconstruction, 
we present a mock data analysis, inspired by the spectra encountered in [11]. A discrete mock 
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Figure 2: (left) Reconstruction of the mock spectral function (gray dashed) using Bryan's prescription on 
N T = 12 datapoints. Note that even the lowest lying peak is washed out and higher lying structures are not 
captured at all. (right) Based on the same data, the MEM in the extended search space 7V e xt = 50 succeeds 
in reconstructing the lowest lying peak. Wiggling behavior at larger co is an artifact of the chosen prior and 
can be identified as such by variation of m(co) 



spectrum is calculated from an analytic expression containing the sum of four Gaussian peaks as 
shown through the gray dashed lines in Fig. 2 (One exponentially small peak is present at negative 
frequencies CO = —2.6). From it, a set of yV T = 12 data-points is generated via eq.(1.2), which after 
being perturbed by Gaussian noise with <jj = 0.01 jDf ea \ j G [1, • • • ,N Z ] is fed to the MEM code 
as input. For the prior we choose the function m(oo) = m _ , so that m{co m { n ) = 0.01. 

On the left of Fig.2 we show that according to Bryan's SVD basis the spectrum is not well 
reproduced. The position and width of even the lowest lying peak is not captured adequately, 
structure at higher CO is completely absent. Unsurprisingly the probability of this solution is very 
small (<2 = log[PMEM[p\D,m]] ~ — 10 4 ). In contrast, after extending the basis to yV ex t = 50, we are 
able to reconstruct from the same N r = 12 dataset the lowest lying peak to a much better degree 
(<2 = —5.5) with both position and width being closer to the mock spectrum. This result is a 
direct counterexample to the claim that the global extremum of Pmem [p \D, m] will always lie in the 
singular subspace. 

At higher frequencies a lot of wiggly structures appear on the right of Fig.2, which are an 
artifact of the choice of prior function. These can be identified as such by repeating the MEM 
reconstruction with different prior functions m(ffl) and observing that only the lowest lying peak 
remains invariant. 

3. Fourier Basis 

We have seen that the MEM solution, i.e. the global maximum of the functional Pmem [p\D, m] 
does not in general lie in the singular subspace obtained from the first N t columns of the SVD 
matrix U. We use this freedom to propose a completely different basis for the search space, which 
cures the central disadvantage of Bryan's SVD basis functions. Our choice decouples the success 
of the reconstruction from the the position of the spectral features relative to t Pmin 

Based on the real Fourier basis,we can parametrize the spectral function in closed form 



pi = m/exp 



Wexl/2 Wext/2 

bi+ b 2 jsm[(coi-co min )j]+ b 2 j+icos[(coi -tiw)j] 
7=1 7=i 
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Figure 3: (left) N T = 12 mock data coming from a single delta peak positioned at increasing frequency 
do- We choose to perturb the data with Gaussian noise, such that a constant relative error ensues, (right) 
Reconstructed peak width versus peak position. The SVD basis yields acceptable results while the peak 
lies within the region, where oscillations exist in the basis functions. Once in the damping region 0) > 5 
the resolution of the reconstructed width deteriorates. (Note that here co 6 [—5,20] so that the SVD basis 
functions oscillate until o osc ~ 5, different from the ones in Fig.l) In contrast, the Fourier basis provides a 
position independent resolution, which even slightly improves due to the better absolute errors in the data at 
larger (Oq. 



First note that since we have an analytic formula available in this case, there is no need to perform 
the SVD to determine the basis functions. Especially for large numbers of N ext > 50 this translates 
into a numerically more robust determination of the spectrum, since for the SVD to accurately 
reproduce the highly oscillating parts of the basis functions a very large numerical precision is 
usually required. 

The main feature of this parametrization is that the choice of co^in becomes irrelevant to the 
success of the reconstruction. To illustrate this point we perform several mock data analyses based 
on single delta peaks placed at different frequencies (Oq. The data obtained from these spectra 
(see Fig. 3 left) shows a simple exponential falloff and is perturbed by Gaussian noise with a a, = 
0.000lDf eal before being fed to the MEM as input. 

As expected from the finite number of datapoints and the finite error present, the single delta 
peak is not fully resolved (i.e. F/Aco > 1), no matter where it is positioned at positive frequencies. 
Fig. 3 shows the reconstructed width versus the peak position, which allows us to observe that the 
MEM resolution deteriorates in the SVD basis for peaks at the higher frequency positions. If the 
peaks lie within the oscillating region of the SVD basis functions their reconstruction only suffers 
mildly from the peak position, but once we position the peak in the purely damped regime, the 
artificially induced width in the reconstruction increases strongly. 

The Fourier basis behaves much more favorably, as it does not show any decrease in resolution 
with change in the peak position. Note that the width of the reconstruction even becomes slightly 
smaller for peaks at high frequencies, since there the absolute error in the mock data is very small. 
Hence we have found a remedy to the central problem connected with Bryan's basis, i.e. the Fourier 
basis allows us to reconstruct peak structures independently from the choice of ftWn- 
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4. Summary and concluding remarks 

We have shown by an inspection of the basis functions used in Bryan's approach that in general 
the MEM solution does not lie in the N r dimensional search space obtained by an SVD of the 
transpose kernel. Since the uniqueness of a possible solution in the full N a dimensional space 
was proven on general grounds in [8], we proposed as a first step to systematically enlarge Bryan's 
search space by including more and more columns of the SVD matrix U. The reconstructed spectra, 
presented in Fig. 2, which were based on mock data inspired by QCD conelators found in [11], 
represent an example of where our improvement becomes crucial for the success of the MEM. 

In order to improve the MEM further, we introduced the Fourier basis to parametrize the 
spectral function based on trigonometric functions. The availability of an analytic expression in 
this case renders the SVD superfluous. Hence the numerical results are more stable, as the SVD 
requires a large numerical precision to correctly determine the highly oscillating features of its 
basis functions. Deploying the Fourier basis to reconstruct single delta peak mock-spectra, we find 
that the resolution of the result does not depend on the peak position, in contrast to the SVD basis, 
which shows a deteriorating resolution for peaks at increasing frequencies. 

The source code of the MEM program used to carry out all of the mock data analyses pre- 
sented here can be found on-line at www.scicode.org/ExtMEM together with a short manual and 
instructions for compilation. 
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